#########################################
## CONJOINT ANALYSIS - UNITED STATES  ###
#########################################


rm(list=ls(all=TRUE))
detach(US.conj)


### Load libraries
library(arm)
library(stargazer)
library(foreign)
library(sjPlot)
library(ggplot2)
library(ggthemes)
library(ggeffects)
library(mediation)
library(cregg)
library(FindIt)


### Read and attach data

#Find your path

US.conj <- read.csv("C:/Users/gmagni/OneDrive/SideProjects/3-AndyNEW/00-Outlin&Paper-Appendix/00-Latest/00-JOP/R&R/FinalSubmission/Data/US-LGBT-Conj.csv")


attach(US.conj)
n = nrow(US.conj)




### Code variables



## Conjoint variables: candidate choice

table(CJ.1.ChoiceDummy)
dummy.choice.B.us <- as.numeric(as.character(CJ.1.ChoiceDummy))
table(dummy.choice.B.us)

table(CJ.1.REC_1)
race.B.us = rep(NA,n)
race.B.us[CJ.1.REC_1=="Asian"]="Asian"
race.B.us[CJ.1.REC_1=="Black"]="Black"
race.B.us[CJ.1.REC_1=="White"]="White"
race.B.us[CJ.1.REC_1=="Latino"]="Latino"
race.B.us[CJ.1.REC_1=="Native American"]="Native American"
table(race.B.us)

race.B.us.factor <- factor(race.B.us, levels = c("White", "Asian", "Black", "Latino", "Native American"))
table(race.B.us.factor)


table(CJ.1.REC_2)
sexorient.B.us = rep(NA,n)
sexorient.B.us[CJ.1.REC_2=="Straight"]="Straight"
sexorient.B.us[CJ.1.REC_2=="Gay"]="Gay"
table(sexorient.B.us)

sexorient.B.us.factor <- factor(sexorient.B.us, levels = c("Straight", "Gay"))
table(sexorient.B.us.factor)


table(CJ.1.REC_3)
religion.B.us = rep(NA,n)
religion.B.us[CJ.1.REC_3=="Not religious"]="Not religious"
religion.B.us[CJ.1.REC_3=="Christian"]="Christian"
religion.B.us[CJ.1.REC_3=="Muslim"]="Muslim"
religion.B.us[CJ.1.REC_3=="Jewish"]="Jewish"
table(religion.B.us)

religion.B.us.factor <- factor(religion.B.us, levels = c("Not religious", "Christian", 
                                                         "Muslim", "Jewish"))
table(religion.B.us.factor)


table(CJ.1.REC_4)
age.B.us = rep(NA,n)
age.B.us[CJ.1.REC_4=="35"]="35"
age.B.us[CJ.1.REC_4=="44"]="44"
age.B.us[CJ.1.REC_4=="56"]="56"
age.B.us[CJ.1.REC_4=="71"]="71"
table(age.B.us)

age.B.us.factor <- factor(age.B.us, levels = c("35", "44", "56", "71"))
table(age.B.us.factor)


table(CJ.1.REC_5)
health.B.us = rep(NA,n)
health.B.us[CJ.1.REC_5=="In a wheelchair since birth"]="Weelchair since birth"
health.B.us[CJ.1.REC_5=="Healthy"]="Healthy"
health.B.us[CJ.1.REC_5=="Overweight, has diabetes"]="Overweight, diabetes"
health.B.us[CJ.1.REC_5=="HIV positive"]="HIV+"
health.B.us[CJ.1.REC_5=="HIV positive since birth"]="HIV+ since birth"
table(health.B.us)

health.B.us.factor <- factor(health.B.us, levels = c("Healthy", "HIV+", "Weelchair since birth", 
                                                     "Overweight, diabetes",
                                                     "HIV+ since birth"))
table(health.B.us.factor)


table(CJ.1.REC_6)
polexp.B.us = rep(NA,n)
polexp.B.us[CJ.1.REC_6=="No previous experience"]="No experience"
polexp.B.us[CJ.1.REC_6=="Member of state legislature"]="State legislature"
polexp.B.us[CJ.1.REC_6=="Member of the U.S. House of Representatives"]="US House of Repres"
table(polexp.B.us)

polexp.B.us.factor <- factor(polexp.B.us, levels = c("No experience", 
                                                     "State legislature", 
                                                     "US House of Repres"))
table(polexp.B.us.factor)


table(CJ.1.REC_7)
educ.B.us = rep(NA,n)
educ.B.us[CJ.1.REC_7=="Less than high school"]="No high school"
educ.B.us[CJ.1.REC_7=="High school degree"]="High school"
educ.B.us[CJ.1.REC_7=="College degree"]="Bachelor"
educ.B.us[CJ.1.REC_7=="Master degree"]="Master"
table(educ.B.us)

educ.B.us.factor <- factor(educ.B.us, levels = c("No high school", "High school", 
                                                 "Bachelor", "Master"))
table(educ.B.us.factor)


table(CJ.1.REC_8)
gender.B.us = rep(NA,n)
gender.B.us[CJ.1.REC_8=="Man"]="Man"
gender.B.us[CJ.1.REC_8=="Woman"]="Woman"
gender.B.us[CJ.1.REC_8=="Transgender"]="Transgender"
table(gender.B.us)

gender.B.us.factor <- factor(gender.B.us, levels = c("Man", "Transgender", "Woman"))
table(gender.B.us.factor)


## Causal mechanisms 

table(CJ.1.CM_CJ.B1.CM_1.Dummy) #left-leaning
cm.left.leaning = rep(NA,n)
cm.left.leaning[CJ.1.CM_CJ.B1.CM_1.Dummy=="0"]=0
cm.left.leaning[CJ.1.CM_CJ.B1.CM_1.Dummy=="1"]=1
table(cm.left.leaning)

table(CJ.1.CM_CJ.CM_5.Dummy) #prefer as neighbor
cm.pref.neighb = rep(NA,n)
cm.pref.neighb[CJ.1.CM_CJ.CM_5.Dummy=="0"]=0
cm.pref.neighb[CJ.1.CM_CJ.CM_5.Dummy=="1"]=1
table(cm.pref.neighb)

table(CJ.1.CM_CJ.CM_6.Dummy) #better chances to win election
cm.electorab = rep(NA,n)
cm.electorab[CJ.1.CM_CJ.CM_6.Dummy=="0"]=0
cm.electorab[CJ.1.CM_CJ.CM_6.Dummy=="1"]=1
table(cm.electorab)



## Controls


table(IncomeH)
income.h.us = rep(NA,n)
income.h.us[IncomeH=="No income"]=0
income.h.us[IncomeH=="Less than  $5.000"]=1
income.h.us[IncomeH=="$5.000-$10.000"]=2
income.h.us[IncomeH=="$10.000-$15.000"]=3
income.h.us[IncomeH=="$15.000-$20.000"]=4
income.h.us[IncomeH=="$20.000-$25.000"]=5
income.h.us[IncomeH=="$25.000-$30.000"]=6
income.h.us[IncomeH=="$30.000-$40.000"]=7
income.h.us[IncomeH=="$40.000-$50.000"]=8
income.h.us[IncomeH=="$50.000-$60.000"]=9
income.h.us[IncomeH=="$60.000-$70.000"]=10
income.h.us[IncomeH=="$70.000-$85.000"]=11
income.h.us[IncomeH=="$85.000-$100.000"]=12
income.h.us[IncomeH=="$100.000-$150.000"]=13
income.h.us[IncomeH=="$150.000-$200.000"]=14
income.h.us[IncomeH=="$200.000-$250.000"]=15
income.h.us[IncomeH=="More than $250.000"]=16
table(income.h.us)

table(Educ)
educ.us = rep(NA,n)
educ.us[Educ=="No high school"]=0
educ.us[Educ=="Some high school but no diploma"]=1
educ.us[Educ=="High school graduate"]=2
educ.us[Educ=="Some college but no degree"]=3
educ.us[Educ=="Associate degree in college"]=4
educ.us[Educ=="Bachelor's degree"]=5
educ.us[Educ=="Master's degree, professional school degree (for example: md, jd), or PhD"]=6
table(educ.us)

table(GenderID)
female.us = rep(NA,n)
female.us[GenderID=="Woman"]=1
female.us[GenderID=="Man"]=0
table(female.us)
prop.table(table(female.us))

table(YOB)
age.us<-2018-YOB
table(age.us)
mean(age.us, na.rm=TRUE)
median(age.us, na.rm=TRUE)


table(PartyID)
party.ID = rep(NA,n)
party.ID[PartyID=="Democrat"]=1
party.ID[PartyID=="Independent"]=2
party.ID[PartyID=="Republican"]=3
table(party.ID)

party.ID.names = rep(NA,n)
party.ID.names[PartyID=="Democrat"]="Dem"
party.ID.names[PartyID=="Independent"]="Rep"
party.ID.names[PartyID=="Republican"]="Ind"
table(party.ID.names)


table(VoteChoice)
party.choice.us = rep(NA,n)
party.choice.us[VoteChoice=="Hillary Clinton"]="Voted Clinton"
party.choice.us[VoteChoice=="Donald Trump"]="Voted Trump"
#party.choice.us[VoteChoice=="Other"]=3
table(party.choice.us)
prop.table(table(party.choice.us))

voted.trump = rep(NA,n) #based only on Trump vs. Clinton voters
voted.trump[party.choice.us=="Voted Clinton"]=0
voted.trump[party.choice.us=="Voted Trump"]=1
table(voted.trump)

table(VoteChoice)
party.choice.2.us = rep(NA,n)
party.choice.2.us[VoteChoice=="Hillary Clinton"]=1
party.choice.2.us[VoteChoice=="Donald Trump"]=2
party.choice.2.us[VoteChoice=="Other"]=3
party.choice.2.us[Vote.Y.N=="No, I didn't vote"]=4
table(party.choice.2.us)
prop.table(table(party.choice.2.us))

voted.trump.all <- as.numeric(party.choice.2.us=="2") #based on all voters
table(voted.trump.all)

table(PolIdeology)
conserv = rep(NA,n)
conserv[PolIdeology=="Extremely liberal"]=1
conserv[PolIdeology=="Liberal"]=2
conserv[PolIdeology=="Slightly liberal"]=3
conserv[PolIdeology=="Moderate or middle of the road"]=4
conserv[PolIdeology=="Slightly conservative"]=5
conserv[PolIdeology=="Conservative"]=6
conserv[PolIdeology=="Extremely conservative"]=7
table(conserv)

table(KnowLGBT)
know.lgbt = rep(NA,n)
know.lgbt[KnowLGBT=="Yes"]=1
know.lgbt[KnowLGBT=="No"]=0
table(know.lgbt)

table(SexOrient)
lgbt.self = rep(NA,n)
lgbt.self[SexOrient=="Gay" | SexOrient=="Lesbian" | SexOrient=="Bisexual" | SexOrient=="Other"]=1
lgbt.self[SexOrient=="Straight"]=0
table(lgbt.self)

table(ReligY.N)
religious.yes = rep(NA,n)
religious.yes[ReligY.N=="Yes"]=1
religious.yes[ReligY.N=="No"]=0
table(religious.yes)

table(ReligID)
relig.ID = rep(NA,n)
relig.ID[ReligID=="Evangelical Protestant"]="Evangelical"
relig.ID[ReligID=="Jewish"]="Jewish"
relig.ID[ReligID=="Mainline Protestant"]="Main Protestant"
relig.ID[ReligID=="Orthodox Christian"]="Orthodox"
relig.ID[ReligID=="Roman Catholic"]="Catholic"
relig.ID[ReligID=="Muslim"]="Muslim"
relig.ID[ReligID=="Mormon"]="Mormon"
relig.ID[ReligID=="Other"]="Other"
table(relig.ID)

table(Religiosity)
religiosity = rep(NA,n)
religiosity[Religiosity=="Never"]=0
religiosity[Religiosity=="Only on special holidays"]=1
religiosity[Religiosity=="At least once a month"]=2
religiosity[Religiosity=="Once a week"]=3
religiosity[Religiosity=="Every day"]=4
table(religiosity)


table(Race)
prop.table(table(Race))
race = rep(NA,n)
race[Race=="White"]="white"
race[Race=="Black or African American"]="black"
race[Race=="Black or African American,White"]="black"
race[Race=="Hispanic/Latino"]="latino"
race[Race=="Hispanic/Latino,White"]="latino"
race[Race=="Asian"]="asian"
race[Race=="Asian,White"]="asian"
race[Race=="American Indian or Alaska Native"]="native"
race[Race=="Native Hawaiian or Pacific Islander"]="native"
race[Race=="American Indian or Alaska Native,White"]="native"
race[Race=="Native Hawaiian or Pacific Islander,White"]="native"
table(race)
prop.table(table(race))

race.factor <- factor(race, levels = c("white", "black", "latino", "asian", "native"))
table(race.factor)


## Add coded variables to dataset

US.conj$soc.conserv.us <- soc.conserv.us
US.conj$income.h.us <- income.h.us
US.conj$educ.us <- educ.us
US.conj$female.us <- female.us
US.conj$age.us <- age.us
US.conj$party.choice.us <- party.choice.us
US.conj$party.choice.2.us <- party.choice.2.us
US.conj$att.check.passed <- att.check.passed
US.conj$party.ID <- party.ID
US.conj$relig.ID <- relig.ID
US.conj$conserv <- conserv
US.conj$race <- race
US.conj$voted.trump <- voted.trump


US.conj$rating.B.us <- rating.B.us
US.conj$dummy.choice.B.us <- dummy.choice.B.us
US.conj$age.B.us <- age.B.us
US.conj$age.B.us.factor <- age.B.us.factor
US.conj$religion.B.us <- religion.B.us
US.conj$religion.B.us.factor <- religion.B.us.factor
US.conj$polexp.B.us <- polexp.B.us
US.conj$polexp.B.us.factor <- polexp.B.us.factor
US.conj$health.B.us <- health.B.us
US.conj$health.B.us.factor <- health.B.us.factor
US.conj$sexorient.B.us <- sexorient.B.us
US.conj$sexorient.B.us.factor <- sexorient.B.us.factor
US.conj$educ.B.us <- educ.B.us
US.conj$educ.B.us.factor <- educ.B.us.factor
US.conj$race.B.us <- race.B.us
US.conj$race.B.us.factor <- race.B.us.factor
US.conj$gender.B.us <- gender.B.us
US.conj$gender.B.us.factor <- gender.B.us.factor
US.conj$know.lgbt <- know.lgbt
US.conj$lgbt.self <- lgbt.self
US.conj$religious.yes <- religious.yes
US.conj$religiosity <- religiosity
US.conj$cm.left.leaning <- cm.left.leaning
US.conj$cm.pref.neighb <- cm.pref.neighb
US.conj$cm.electorab <- cm.electorab




#################

### Analysis




# General model

conjoint.us.lgbt <- lm(dummy.choice.B.us ~ gender.B.us.factor + sexorient.B.us.factor + 
                          health.B.us.factor +
                          age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                          race.B.us.factor + educ.B.us.factor, data = US.conj)
summary(conjoint.us.lgbt)


########



### Analysis of LGBT conjoint using [cregg] ###

#Rename conjoint attribute variables to have clearer plots and attach them to dataset 

Sexual_Orientation <- sexorient.B.us.factor
Gender_Identity <- gender.B.us.factor
Education <- educ.B.us.factor
Race <- race.B.us.factor
Political_Experience <- polexp.B.us.factor
Religion <- religion.B.us.factor
Age <- age.B.us.factor
Health <- health.B.us.factor

US.conj$Sexual_Orientation <- Sexual_Orientation
US.conj$Gender_Identity <- Gender_Identity
US.conj$Education <- Education
US.conj$Race <- Race
US.conj$Political_Experience <- Political_Experience
US.conj$Religion <- Religion
US.conj$Age <- Age
US.conj$Health <- Health

US.conj$cm.electorab <- cm.electorab
US.conj$cm.left.leaning <- cm.left.leaning
US.conj$cm.pref.neighb <- cm.pref.neighb


Choice <- dummy.choice.B.us
US.conj$Choice <- Choice


#Create dummy variables for subgroup analysis and attach them to dataset

Democrats = rep(NA,n)
Democrats[party.ID==1]=1
Democrats[party.ID==3]=0 #conserv
table(Democrats)
US.conj$Democrats <- Democrats

Liberal = rep(NA,n)
Liberal[conserv==1 | conserv==2]=1
Liberal[conserv==6 | conserv==7]=0
table(Liberal)
US.conj$Liberal <- Liberal

LGBT.friends <-know.lgbt
table(LGBT.friends)
US.conj$LGBT.friends <- LGBT.friends

LGBT.voters <- lgbt.self
table(lgbt.self)
US.conj$lgbt.self <- lgbt.self


Women<-female.us
table(Women)
US.conj$Women <- Women

Young = rep(NA,n)
Young[age.us<="35"]=1
Young[age.us>="60"]=0
table(Young)
US.conj$Young <- Young

Non.religious = rep(NA,n)
Non.religious[religiosity==0]=1
Non.religious[religiosity==3 | religiosity==4]=0
table(Non.religious)
US.conj$Non.religious <- Non.religious

table(income.h.us)
Low.income = rep(NA,n)
Low.income[income.h.us==0 | income.h.us ==2 | income.h.us==3 | income.h.us==4 | income.h.us==5]=1
Low.income[income.h.us!=0 & income.h.us !=2 & income.h.us!=3 & income.h.us!=4 & income.h.us!=5]=0
table(Low.income)
US.conj$Low.income <- Low.income

High.income = rep(NA,n)
High.income[income.h.us==13 | income.h.us ==14 | income.h.us==15 | income.h.us==16]=1
High.income[income.h.us!=13 & income.h.us !=14 & income.h.us!=15 & income.h.us!=16]=0
table(High.income)
US.conj$High.income <- High.income

table(race)
Whites = rep(NA,n)
Whites[race=="white"]=1
Whites[race!="white"]=0
table(Whites)
US.conj$Whites <- Whites

Blacks = rep(NA,n)
Blacks[race=="black"]=1
Blacks[race!="black"]=0
table(Blacks)
US.conj$Blacks <- Blacks

Latinx = rep(NA,n)
Latinx[race=="latino"]=1
Latinx[race!="latino"]=0
table(Latinx)
US.conj$Latinx <- Latinx

Asians = rep(NA,n)
Asians[race=="asian"]=1
Asians[race!="asian"]=0
table(Asians)
US.conj$Asians <- Asians


## Subset of dataset 

US.conj.republ <- subset(US.conj, party.ID == 3)
US.conj.dem <- subset(US.conj, party.ID == 1)
US.conj.conserv <- subset(US.conj, conserv == 6 | conserv == 7)
US.conj.liber <- subset(US.conj, conserv == 1 | conserv == 2)
US.conj.trump <- subset(US.conj, party.choice.us == "Voted Trump")
US.conj.clinton <- subset(US.conj, party.choice.us == "Voted Clinton")
US.conj.evangel <- subset(US.conj, relig.ID == "Evangelical")
US.conj.women <- subset(US.conj, female.us == 1)
US.conj.men <- subset(US.conj, female.us == 0)
US.conj.knowlgbt <- subset(US.conj, know.lgbt == 1)
US.conj.notknowlgbt <- subset(US.conj, know.lgbt == 0)
US.conj.noreligserv <- subset(US.conj, religiosity == 0)
US.conj.religserv <- subset(US.conj, religiosity == 3 | religiosity == 4)
US.conj.young <- subset(US.conj, age.us <= 35)
US.conj.old <- subset(US.conj, age.us >= 60)
US.conj.white <- subset(US.conj, race == "white")
US.conj.black <- subset(US.conj, race == "black")
US.conj.latino <- subset(US.conj, race == "latino")
US.conj.asian <- subset(US.conj, race == "asian")
US.conj.native <- subset(US.conj, race == "native")
US.conj.job.sec <- subset(US.conj, job.insec.us == 0)
US.conj.job.insec <- subset(US.conj, job.insec.us == 2 | job.insec.us == 3 | job.insec.us == 4)
US.conj.lgbtvoters <- subset(US.conj, lgbt.self == 1)
US.conj.straightvoters <- subset(US.conj, lgbt.self == 0)
US.conj.lowincome <- subset(US.conj, income.h.us==0 | income.h.us ==2 | income.h.us==3 | income.h.us==4 | income.h.us==5)
US.conj.highincome <- subset(US.conj, income.h.us==13 | income.h.us ==14 | income.h.us==15 | income.h.us==16)


US.conj.liber.LGBTfriends <- subset(US.conj, know.lgbt == 1 & (conserv == 1 | conserv == 2))
US.conj.cons.LGBTfriends <- subset(US.conj, know.lgbt == 1 & (conserv == 6 | conserv == 7))


table(lgbt.self)


### ACME ###

#ACME - General
amce(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
       Health + Age + Religion + Political_Experience + 
       Race + Education, id = ~ ResponseId)

amce.1 <- amce(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)


plot(amce.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())




### MM ###

mm.1 <- mm(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
     Health + Age + Religion + Political_Experience + 
     Race + Education, id = ~ ResponseId, h0 = 0.5)

plot(mm.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) + geom_vline(xintercept = 0.50, colour="gray33")



## SUBGROUP ANALYSIS

# AMCE by groups (just change dataset in command)

amce.sub <- amce(US.conj.men, Choice ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

amce.sub <- amce(US.conj.women, Choice ~ Sexual_Orientation + Gender_Identity +  
                          Health + Age + Religion + Political_Experience + 
                          Race + Education, id = ~ ResponseId)

# MMs split by groups

#MM by gender
stacked.gender <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ Women)
plot(stacked.gender, group = "Women", feature_headers = FALSE) + geom_vline(xintercept = 0.50)

#MM by party
stacked.party <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                      Health + Age + Religion + Political_Experience + 
                      Race + Education, id = ~ ResponseId,
                    estimate = "mm", by = ~ Democrats)

#MM by political ideology
stacked.polideology <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, id = ~ ResponseId,
                          estimate = "mm", by = ~ Liberal)

#MM by having LGBT friends
stacked.lgbtfriends <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, id = ~ ResponseId,
                          estimate = "mm", by = ~ LGBT.friends)

#MM by age
stacked.age <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                    Health + Age + Religion + Political_Experience + 
                    Race + Education, id = ~ ResponseId,
                  estimate = "mm", by = ~ Young)

#MM by religiosity
stacked.religiosity <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, id = ~ ResponseId,
                          estimate = "mm", by = ~ Non.religious)

#MM by sexual orientation of voters
stacked.sexor <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, id = ~ ResponseId,
                          estimate = "mm", by = ~ lgbt.self)

#MM by race
stacked.race <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                      Health + Age + Religion + Political_Experience + 
                      Race + Education, id = ~ ResponseId,
                    estimate = "mm", by = ~ race)

#MM by low income (combine it with high income below)
stacked.lowincome <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                      Health + Age + Religion + Political_Experience + 
                      Race + Education, id = ~ ResponseId,
                    estimate = "mm", by = ~ Low.income)

#MM by high income (combine it with low icnome above)
stacked.highincome <- cj(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                      Health + Age + Religion + Political_Experience + 
                      Race + Education, id = ~ ResponseId,
                    estimate = "mm", by = ~ High.income)

### SUBGROUP DIFFERENCES ###

### ACME differences
acme.diff.1 <- amce_diffs(US.conj, dummy.choice.B.us ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, ~Women , id = ~ ResponseId)


plot(acme.diff.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())


### MM differences
  
#mm differences by gender
mm.diff.gender <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                             Health + Age + Religion + Political_Experience + 
                             Race + Education, 
         ~ Women, id = ~ ResponseId)
plot(mm.diff.gender, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by party
mm.diff.party <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                             Health + Age + Religion + Political_Experience + 
                             Race + Education, 
                           ~ Democrats, id = ~ ResponseId)
plot(mm.diff.party, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by political ideology
mm.diff.polideol <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, 
                          ~ Liberal, id = ~ ResponseId)
plot(mm.diff.polideol, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by having LGBT friends
mm.diff.lgbtfriends <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                               Health + Age + Religion + Political_Experience + 
                               Race + Education, 
                             ~ LGBT.friends, id = ~ ResponseId)
plot(mm.diff.lgbtfriends, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by age
mm.diff.age <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                                  Health + Age + Religion + Political_Experience + 
                                  Race + Education, 
                                ~ Young, id = ~ ResponseId)
plot(mm.diff.age, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by religiosity
mm.diff.religiosity <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                          Health + Age + Religion + Political_Experience + 
                          Race + Education, 
                        ~ Non.religious, id = ~ ResponseId)
plot(mm.diff.religiosity, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by sexual orientation
mm.diff.sexor <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                                  Health + Age + Religion + Political_Experience + 
                                  Race + Education, 
                                ~ lgbt.self, id = ~ ResponseId)


#mm differences by race
mm.diff.race <- mm_diffs(US.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, 
                          ~ race, id = ~ ResponseId)





## CAUSAL MECHANISM

# (substitute dataset for subgroup analysis)

# Ideology - left
conjoint.cm.left <- lm(cm.left.leaning ~ gender.B.us.factor + sexorient.B.us.factor + 
                          health.B.us.factor +
                          age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                          race.B.us.factor +
                          educ.B.us.factor, data = US.conj)
summary(conjoint.cm.left)


# Prejudice
conjoint.cm.neighb <- lm(cm.pref.neighb ~ gender.B.us.factor + sexorient.B.us.factor + 
                            health.B.us.factor +
                            age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                            race.B.us.factor +
                            educ.B.us.factor, data = US.conj)
summary(conjoint.cm.neighb)

# Electability
conjoint.cm.elector <- lm(cm.electorab ~ gender.B.us.factor + sexorient.B.us.factor + 
                             health.B.us.factor +
                             age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                             race.B.us.factor +
                             educ.B.us.factor, data = US.conj)
summary(conjoint.cm.elector)

#5B - electab controlling for prejudice
conjoint.cm.elector.prej <- lm(cm.electorab ~ gender.B.us.factor + sexorient.B.us.factor + 
                            health.B.us.factor +
                            age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                            race.B.us.factor + educ.B.us.factor +
                             cm.pref.neighb, data = US.conj)
summary(conjoint.cm.elector.prej)



#Causal mediation

library(mediation)

#New with pre-treatment covariates


med.fit <- lm(cm.pref.neighb ~ gender.B.us.factor + sexorient.B.us.factor + 
                 health.B.us.factor + age.B.us.factor + religion.B.us.factor + 
                 polexp.B.us.factor + race.B.us.factor + educ.B.us.factor +
                 age.us + educ.us + female.us + know.lgbt + party.ID + race +
                 religiosity + income.h.us, data = US.conj)
out.fit <- lm(dummy.choice.B.us ~ cm.pref.neighb + 
                gender.B.us.factor + sexorient.B.us.factor + 
                 health.B.us.factor + age.B.us.factor + religion.B.us.factor + 
                 polexp.B.us.factor + race.B.us.factor + educ.B.us.factor +
                 age.us + educ.us + female.us + know.lgbt + party.ID + race +
                 religiosity + income.h.us, data = US.conj) 

##
med.out <- mediate(med.fit, out.fit, treat = "gender.B.us.factor", 
                   mediator = "cm.pref.neighb", 
                   sims = 1000,
                   dropobs = TRUE,
                   robustSE = TRUE)
summary(med.out)
plot(med.out)


#Sensitivity analysis
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)
plot(sens.out, sens.par = "rho", main = "Meritocracy", ylim = c(-0.2, 0.2))


# re-code for sensitivity analysis for treatment with more than 2 conditions

table(gender.B.us.factor)
treat.sens = rep(NA,n)
treat.sens[gender.B.us.factor=="Man"]=0
treat.sens[gender.B.us.factor=="Transgender"]=1
treat.sens[gender.B.us.factor=="Woman"]=2
table(treat.sens)
treat.sens<-factor(treat.sens)

table(health.B.us.factor)
treat.HIV = rep(NA,n)
treat.HIV[health.B.us.factor=="Healthy"]=0
treat.HIV[health.B.us.factor=="HIV+"]=1
treat.HIV[health.B.us.factor=="Weelchair since birth"]=2
treat.HIV[health.B.us.factor=="Overweight, diabetes"]=3
treat.HIV[health.B.us.factor=="HIV+ since birth"]=4
table(treat.HIV)
treat.HIV<-factor(treat.HIV)


med.fit <- lm(cm.pref.neighb ~ gender.B.us.factor + sexorient.B.us.factor + 
                treat.HIV + age.B.us.factor + religion.B.us.factor + 
                polexp.B.us.factor + race.B.us.factor + educ.B.us.factor +
                age.us + educ.us + female.us + know.lgbt + party.ID + race +
                religiosity + income.h.us, data = US.conj)
out.fit <- lm(dummy.choice.B.us ~ cm.pref.neighb + 
                gender.B.us.factor + sexorient.B.us.factor + 
                treat.HIV + age.B.us.factor + religion.B.us.factor + 
                polexp.B.us.factor + race.B.us.factor + educ.B.us.factor +
                age.us + educ.us + female.us + know.lgbt + party.ID + race +
                religiosity + income.h.us, data = US.conj) 

##
med.out <- mediate(med.fit, out.fit, treat = "treat.HIV", 
                   mediator = "cm.pref.neighb", 
                   sims = 1000,
                   dropobs = TRUE,
                   robustSE = TRUE)
summary(med.out)

#Sensitivity analysis
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)



## INTERACTION EFFECTS: FINDIT PACKAGE

#Subset so you don't have missing data
myvars2 <- c("dummy.choice.B.us", "religion.B.us.factor", "gender.B.us.factor", "sexorient.B.us.factor",
            "health.B.us.factor", "educ.B.us.factor", "age.B.us.factor", "religion.B.us.factor",
            "polexp.B.us.factor", "race.B.us.factor", 
            "ResponseId")
US.conj.LGBT.inter <- US.conj.liber[myvars2]
US.conj.LGBT.inter <- na.omit(US.conj.LGBT.inter)

## Specify the order of each factor
US.conj.LGBT.inter$gender.B.us.factor<- factor(US.conj.LGBT.inter$gender.B.us.factor,ordered=TRUE,
                                        levels=c("Woman", "Transgender", "Man"))
US.conj.LGBT.inter$sexorient.B.us.factor<- factor(US.conj.LGBT.inter$sexorient.B.us.factor,ordered=TRUE,
                                               levels=c("Straight", "Gay"))
US.conj.LGBT.inter$race.B.us.factor<- factor(US.conj.LGBT.inter$race.B.us.factor,ordered=TRUE,
                                                  levels=c("Black", "Asian", "Latino", 
                                                           "Native American", "White"))

#Calculate AMIE 

#Gender and sexual orientation 
fit.1 <- CausalANOVA(formula=dummy.choice.B.us ~ gender.B.us.factor + sexorient.B.us.factor + 
                       health.B.us.factor +
                       age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                       race.B.us.factor +
                       educ.B.us.factor,
                    int2.formula = ~ gender.B.us.factor:sexorient.B.us.factor,
                    data=US.conj.LGBT.inter, 
                    cluster=US.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.1)


#Sexual orientation and race
fit.2 <- CausalANOVA(formula=dummy.choice.B.us ~ gender.B.us.factor + sexorient.B.us.factor + 
                       health.B.us.factor +
                       age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                       race.B.us.factor +
                       educ.B.us.factor,
                     int2.formula = ~ sexorient.B.us.factor:race.B.us.factor,
                     data=US.conj.LGBT.inter, 
                     cluster=US.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.2)


#Gender and race
US.conj.LGBT.inter$gender.B.us.factor<- factor(US.conj.LGBT.inter$gender.B.us.factor,ordered=TRUE,
                                               levels=c("Woman", "Man", "Transgender"))
fit.3 <- CausalANOVA(formula=dummy.choice.B.us ~ gender.B.us.factor + sexorient.B.us.factor + 
                       health.B.us.factor +
                       age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                       race.B.us.factor +
                       educ.B.us.factor,
                     int2.formula = ~ gender.B.us.factor:race.B.us.factor,
                     data=US.conj.LGBT.inter, 
                     cluster=US.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.3)


#Health and race
US.conj.LGBT.inter$health.B.us.factor<- factor(US.conj.LGBT.inter$health.B.us.factor,ordered=TRUE,
                                               levels=c("HIV+ since birth", "Weelchair since birth", 
                                                        "Healthy", "Overweight, diabetes", "HIV+"))
fit.4 <- CausalANOVA(formula=dummy.choice.B.us ~ gender.B.us.factor + sexorient.B.us.factor + 
                       health.B.us.factor +
                       age.B.us.factor + religion.B.us.factor + polexp.B.us.factor + 
                       race.B.us.factor +
                       educ.B.us.factor,
                     int2.formula = ~ health.B.us.factor:race.B.us.factor,
                     data=US.conj.LGBT.inter, 
                     cluster=US.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.4)



